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ABSTRACT 

In an attempt to understand the extraordinarily small mass-loss rates of late-type O dwarfs, mass fluxes in the relevant 
part of {T e ff, gr)-space are derived from first principles using a previously-described code for constructing moving 
reversing layers. From these mass fluxes, a weak-wind domain is identified within which a star's rate of mass loss 
by a radiatively-driven wind is less than that due to nuclear burning. The five weak-wind stars recently analysed by 
Marcolino et al. (2009) fall within or at the edge of this domain. But although the theoretical mass fluxes for these 
stars are w 1.4 dex lower than those derived with the formula of Vink et al. (2000), the observed rates are still not 
matched, a failure that may reflect our poor understanding of low-density supersonic outflows. 

Mass fluxes are also computed for two strong- wind 04 stars analysed by Bouret et al. (2005). The predictions agree 
with the sharply reduced mass loss rates found when Bouret et al. take wind clumping into account. 
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1. Introduction 

Investigations of late-type O dwarfs in the SMC (Bouret et 
al. 2003; Martins et al. 2004) and the Galaxy (Marcolino et 
al. 2009; M09) find that mass-loss rates $ determined by 
spectral modelling are at least an order of magnitude lower 
than those predicted with the formula of Vink et al. 

(2000). 

Although this weak wind problem is considered to be 
a major problem for stellar wind theory (Puis, Vink, & 
Najarro 2008), this may not in fact be so since the <Ev's 
are semi-empirical estimates based on Monte Carlo (MC) 
calulations of the deposition rates of radiative energy in the 
supersonic winds. Accordingly, this conflict might be a con- 
sequence of the Vink et al. assumptions concerning physical 
conditions in these winds. Moreover, the global dynamical 
constraint imposed by Vink et al and previously by Abbott 
& Lucy (1985) does not guarantee that the derived <&'s are 
consistent with stationary transonic flows. 

In order, therefore, to make a more decisive confronta- 
tion between theory and observation, this paper reports an 
extensive grid of mass fluxes J = &/4irR 2 for O stars us- 
ing a code described earlier (Lucy 2007; L07a). This code, 
which updates the moving reversing layer (RL) theory of 
Lucy & Solomon (1970; LS70), uses a MC technique to es- 
timate the J that allows the flow to accelerate continuously 
from sub- to supersonic velocities. In other words, the eigen- 
value J is determined from first principles by imposing a 
regularity condition at the sonic point. 

2. Solution technique 

Because the code is unchanged from L07a, details are omit- 
ted. However, the search procedure for the adjustable pa- 
rameters has been revised to facilitate a survey of (T e ff,g)- 
space. 

Send offprint requests to: L.B.Lucy 



2.1. Implicit assumptions 

An unstated assumption in L07a is that the subsonic layers 
of a RL are affected by the exterior supersonic layers only 
through backscattered radiation. Specifically, therefore, no 
structural changes occur due to information propagating 
back into the RL via radiative-acoustic waves (Abbott 
1980). This assumption was checked in L07b and found to 
be justified: for the models from L07a matching supersonic 
solutions were constructed and the group velocity of Abbott 
waves computed. In each case, the direction of information 
propagation is everywhere outwards towards higher veloci- 
ties. 

However, even without this justification, there are two 
further reasons for adopting this assumption. First, any 
conflict with observations would then be possible evidence 
that structural changes do indeed occur. Second, since 
such waves have to propagate back into the photosphere 
through an outflow that, according to most spectroscopists, 
is chaotic and clumpy, their information content might well 
be scrambled and lost. 

A second unstated assumption is that after achieving 
supersonic velocities matter continues to be accelerated at 
least until local escape velocity is reached, so that no matter 
falls back onto the RL. Again, for the models of L07a, this 
is confirmed by the matching exterior solutions reported in 
L07b. But this assumption might well be violated for some 
late- type O stars (Howk et al. 2000). 

In addition, note that these dynamic RL's are con- 
structed without the Sobolev approximation, which, for a 
microturbulent velocity of lOfcms^ 1 , is not valid for Mach 
numbers Si 3 - see Fig. 4 in L07a - and so is inappropriate 
for the sub- and transonic flows of O stars. 
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2.2. Adjustable parameters 

In addition to J, a model RL as described in L07a has two 
further adjustable parameters S and s. These define the 
velocity dependence of ge, the radiative acceleration due to 
lines, according to the formula 



ge = g*mcix 



(^j with g* = g-g e 



(1) 



where g e is the radiative acceleration due to electron scat- 
tering. With < 5 < 1 and s > 0, this formula is such that 
the effective gravity g e // = .9 — 5e — ge — at v = a and 
is negative at supersonic velocities, thus allowing the rapid 
acceleration of the transonic flow to be accurately modelled 
- see Fig. 2 in L07a. 

2.3. An identity 

If the equation of motion for plane-parallel, radiatively- 
driven isothermal flow - Eq. (2) in L07a - is integrated 
between heights x\ and a; 2 , we find that <2i j2 = 0, where 



1 . o q. . TT12 , 1 f 2 

2 mi a 2 y xi 



(2) 



Here m = u/a is the flow's Mach number relative to the 
isothermal speed of sound a. 

In physical terms, Q1.2 = implies that the gain in me- 
chanical energy between x\ and X2 is due to the work done 
by the gradients of gas and radiation pressures. But for 
a solution generated with the L07a code, the identity will 
in general be violated because of a non-optimum parame- 
ter vector (J, S, s). Nevertheless, by adjusting this vector as 
described below, we can find the solution giving Q iy2 = 0. 

The above identity holds for any pair of points (xi, x 2 ). 
Here, because of our interest in finding a solution repre- 
senting smooth transonic flow, we choose these points to 
be where mi « 0.5 and m 2 « 2.0, thus straddling the sonic 
point. Accordingly, for this MC technique, the constraint 
Qi,2 = for the interval (mi,m 2 ) is the analogue of the 
regularity condition at v — a in conventional investigations 
of transonic flow. 



2.4. Determining the eigenvalue J 

The steps followed in obtaining a satisfactory model are: 

1) For the chosen parameters T e ff,g, continuum fluxes 
are extracted from the corresponding TLUSTY model 
(Lanz & Hubeny 2003) at 26 frequencies Vk{Hz) from 14.5- 
16.2 dex. This is illustrated in Fig. 1, with parameters ap- 
propriate for a late-type O dwarf. In the MC calculation, 
fluxes at the base of the RL are obtained by linear loga- 
rithmic interpolation in the intervals (v^, Vk+i)- Note that 
because the emergent MC flux at the top of the RL is con- 
strained by Eq. (7) of L07a to = crT^j, only relative con- 
tinuum fluxes of the TLUSTY models are used. 

2) For the same TLUSTY model, ground-state depar- 
ture coefficients are extracted from the point where T e w 
0.75T e ff in order to compute ionization in the RL from Eq. 
(6) of L07a. 

3) Initial values are selected for the parameters (J, S, s), 
and the resulting stratification of the RL computed as de- 
scribed in Sect. 2.3 of L07a. 




15.8 15.6 15.4 15.2 

Log V (llz) 

Fig. 1. Emergent spectrum of TLUSTY model atmosphere with 
parameters as indicated. The open circles locate the points used 
to approximate the continuum. 



4) The radiation field throughout the RL is then derived 
with the MC technique described in Sect. 3.2 of L07a. In 
particular, the MC estimator ge for the radiative accelera- 
tion due to lines is evaluated. 

5) With these values ge , the integral in Eq. (2) is ap- 
proximated as a summation and Qi j2 calculated. 

6) Steps 3) to 6) are repeated for different J's but fixed 
(6, s) in order to locate the intersection with Qi :2 = - 
see Fig. 2. The intersection is derived from a least squares 
linear fit to 5-12 Q- values closely bracketing Q = 0. 

7) From the model thus derived, fits to the variations of 
ge in the sub- and supersonic flows suggest improved values 
of d and s, respectively. If these differ significantly from the 
current values, steps 3) to 7) are repeated. 

2.5. Dynamical consistency 

A qualitative consistency check is the mismatch between 
ge, the computed and ge(v; 5, s), the assumed radiative ac- 
celerations due to lines. These are plotted for the model 
with T e ff = 32, 500K in Fig. 3, which may be compared 
to the plot for T e ff = 50, 000K in L07a. This comparison 
shows that the adopted functional form is here less success- 
ful at capturing the structure of the function ge(v). Clearly, 
a more flexible representation is desirable to achieve better 
dynamical consistency. 

This lack of consistency implies some uncertainty in the 
predicted J's. A rough estimate of this is derived as follows: 
22 of the models in Table 1 were originally computed with 
S = 0.5 and s = 2. When S and s were optimized, the 
changes \AlogJ\ were < 0.05 in 13 cases and > 0.2 in only 
3 cases, with the largest difference being 0.42. 

This experiment shows that the J's are moderately in- 
sensitive to the vector (<5, s) and thus to deviations from 
dynamical consistency. The most important step in deriv- 
ing an accurate J is locating the intersection Q1.2 = 0, 
thus ensuring that radiative driving makes its mandatory 
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Table 1. Computed eigenvalues J(gms 1 cm 2 ) 
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Fig. 2. The search for the eigenvalue J. The atmospheres' pa- 
rameters are T eff = 32500if (filled circles) and T eff = 37500K 
(open circles) both with log g = 4. For each atmosphere, a se- 
quence of models with varying J but fixed (8, s) is plotted. The 
values of log J given by the intersections with Q\.2 = are 
indicated. 
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Fig. 3. Comparison of assumed and computed radiative acceler- 
ations due to lines. The filled circles are the MC estimates and 
the solid line shows the variation assumed in deriving the RL's 
stratification. The open circle is the sonic point. 
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3. Computed mass fluxes 

In this section, J's are calculated for parameter vectors 
(T e ff,g) chosen to illustrate the onset of powerful winds 
for stars. Results are given in Table 1 for 29 models 
with a range of g's appropriate for stars in their H-burning 
main sequence phase. The models' composition is solar with 
N He /N H = 0.1 (Grevess & Sauval 1998) and the included 
metal ions are as in Table 1 of Lanz & Hubeny (2003). 
The microturbulent velocity v t — lOkms^ 1 . In cases where 
dynamical consistency is poor, log J is followed by a colon. 
Note that two models from L07a are recomputed in Table 
1. The changes in logJ are +0.03 for D-50 and -0.07 for 
D-40. 

Throughout this paper, logarithmic values of J are re- 
ported with mass flux unit gm s~ x cmT 2 . If the radius of 
the star is known, its mass-loss rate $ in solar masses per 
year can then be derived from the formula 



R 

log $ = log J + 2 log — h C 



Rr. 



© 



(3) 



with C 



-3.015. 



contribution to the gain of mechanical energy as the flow 
accelerates through the sonic point from mi to 771,2 • 

Given that the disagreement discussed in Sect. 1 is by at 
least 1.0 dex in the uncertainties in J are not of present 
concern. 



3.1. Weak-wind domain 

In thermal equilibrium, a star's mass loss rate via its wind 
equals that by nuclear burning when $ = L/c 2 , where L = 
AttR 2 J-. Accordingly, we can define the weak-wind domain 
in (T e ff , g)-space to be where J < J* = F/c 2 . The locus 
where J = J* has therefore been derived by interpolating 
between models in Table 1 and is shown as the bold line in 
Fig. 4. Also shown is the contour where J/ J* has increased 
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Table 2. Mass fluxes for observed stars 



Fig. 4. Weak-wind domain. The bold line is the locus J = J* 
separating weak from strong winds. Also shown is the locus 
where J exceeds J* by 0.5 dex. The filled circles are the weak- 
wind stars investigated in M09, and the asterisks are ZAMS 
models of Pols et al. (1998) for Z = 0.02 with masses 12(2)34 in 
solar units. 



to 0.5 dex, illustrating the steepness of the onset of strong 
winds. 

The five weak-wind stars analysed in M09 are also plot- 
ted in Fig. 4. They either fall within the weak-wind domain 
or have error bars that just overlap the contour J = J* . 

In order to determine the transition from weak to strong 
winds in terms of stellar mass, zero-age main sequence 
(ZAMS) models from Pols et al. (1998) for Z = 0.02 are 
also included in Fig. 4. From this sequence, we see that stars 
on the ZAMS with M/Mq ^ 22 are in the weak-wind do- 
main while stars with M/Mq <; 28 are in the strong- wind 
domain. 

The corresponding critical masses can be derived for 
the mass- loss formula of Vink et al. (2000). From their Eq. 
(12) with Voo/vesc — 2.6 (Lamers et al. 1995) and noting 
that the formula applies only down to T e ff = 27, 500K, 
we find that the upper limit for weak winds drops from 22 
to < 12, while the lower limit for strong winds drops from 
28 to 14 solar masses. These large changes have profound 
implications for the evolution of massive stars. 

3.2. Individual Marcolino stars 

Finding that the stars of M09 fall in or at the edge of the 
weak-wind domain in Fig. 4 appears to suggest that the J's 
predicted from first principles with the L07a code solves the 
weak wind problem. But before accepting this conclusion, 
we must compare observed and predicted mass fluxes for 
the five stars individually. This is done in Table 2, where 
the data is as follows: 

Column 2: log Jm computed from the measured log $'s 
reported in M09 together with their conservative error es- 
timates. 

Column 3: logJ^ obtained by interpolation in Table 1 
for the T e ff and g given in Table 3 of M09. The computed 
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errors reflect only the propagation of the errors in T e f / and 
g given in M09 and shown in Fig. 4. 

Column 4: logJy derived from Eq. (12) in Vink et al. 
(2000) with parameters from Table 3 in M09. 

Column 5: logj*, where J* = T/c 2 . 

From Table 2, we see that the Jl's are systematically 
smaller than the Jy's by on average 1.4 dex. Nevertheless, 
this huge reduction is still not enough to explain the mea- 
sured values Jm, which are smaller than the J^'s by on 
average 0.8 dex. 

Given the large error bars, the individual differences 
logJu — logJh are barely significant. But with all differ- 
ences being of the same sign, the data as a whole could 
be said to reject the hypothesis that the L07a theory cor- 
rectly predicts observed J's. But this is not purely a ques- 
tion of statistics. The reliability of the diagnostic investiga- 
tions can be questioned, especially for low $'s when often 
only the C IV doublet is available for analysis, and this 
line's weak or absent emission component is unexplained 
(Martins et al. 2004). This latter problem emphasizes how 
limited is our understanding of the supersonic zones of these 
stars' low density winds, a problem affecting the Vink et al. 
(2000) predictions as well as diagnostic codes. 

3.3. The matching problem 

As noted in Scct.l, $'s derived from a global dynamical con- 
straint may be inconsistent with stationary transonic flow; 
and the same remark applies to $'s derived diagnostically 
from circumstellar spectra. To illustrate the magnitude of 
these inconsistencies for the weak- wind stars, the consis- 
tency check carried out in Fig. 3 has been repeated in Fig. 5 
with Jl increased by 1.4 dex to represent the Vink et al 
recipe and decreased by 0.8 dex to represent the Marcolino 
et al. estimates. In both cases, there are huge differences be- 
tween gt(v) and g~e, indicating gross violations of dynamical 
consistency. Moreover, searches in (<5, s)-space with J fixed 
at these high or low J's fail to find alternative dynamically- 
consistent solutions. 

Evidently, any estimate of $ derived from an analysis 
of a wind's supersonic flow should ideally be confirmed by 
matching to a stationary transonic flow, thus exhibiting a 
successful transition to the quasi-static atmosphere where 
the star's absorption line spectrum is formed. 



3.4. Two strong-wind stars 

Given the large discrepancies for the weak-wind stars be- 
tween Jm, Jl and Jy, a reader might reasonably conclude 
that, since no two methods even approximately agree, there 
is no basis for preferring one estimate over another. The 
data in Table 1 is therefore now used to compute J's for the 
two Galactic 04 stars investigated by Bouret et al. (2005). 
These authors are co-authors of the Marcolino paper, and 
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Fig. 5. Departures from dynamical consistency. The solid line 
shows the variation of ge(v) assumed in deriving the strati- 
fication of the RL with the indicated parameters, for which 
Jh{gm cmT 2 s -1 ) = -7.32 dex - see Fig.3. The filled circles are 
the MC estimates ge when J is increased to -5.92 dex; and the 
small open circles are the estimates for J = -8.12 dex. The large 
open circle is the sonic point. 

the same diagnostic codes, TLUSTY and CMFGEN, are 
used. 

For the 04 V star HD 96715, Bouret et al derive $ = 
2.5±0.5 x 10~ 7 Mqut^ 1 , corresponding to logj = — 5.74± 
0.08. This is their preferred estimate, derived with the filling 
factor of clumps as an additional adjustable parameter. For 
a homogeneous wind, $ is 0.86 dex higher, corresponding 
to log J = -4.88. 

The atmospheric parameters derived by Bouret et al. are 
T eff = 43500 ± 2200if and logg = 4.0 ± 0.1. Interpolation 
in Table 1 gives logJh = -5.54 ±0.25. But note that Table 
1 is computed with vt = lOfcms -1 , whereas Bouret et al 
derive v t = 15 ± 5fcms _1 . According to Table 1 in L07a, 
this increase in v t reduces Jl by about 0.30 dex, bringing 
the prediction to logj^ = — 5.84± 0.36. This agrees within 
errors with the clump- corrected value of Bouret et al. but 
is in strong disagreement with the value for a homogeneous 
wind. 

The second Bouret et al. star is the 04 I supergiant 
HD 190429A. The clumped- wind analysis gives $ = 1.8 ± 
0.4 x lO _6 .M02/r _1 , corresponding to logj = —5.31 ± 
0.08. The atmospheric parameters derived by Bouret et al. 
are T eff = 39000 ± 2000K and logg = 3.6 ± 0.1 giving 
logJh = —5.04 ± 0.76, where the large standard error is 
due to the steep gradient of Jl at this star's location in 
(T e //,#)-space. As before, we now apply a correction of - 
0.30 dex because the measured v t = 15±2fcms _1 . The final 
prediction is therefore log Jl = — 5. 34±0. 77, consistent with 
the clumped wind analysis but, in this case, also with the 
value logj = —4.78 ± 0.08 for a homogeneous wind. 

The clumped- wind $'s derived by Bouret et al are lower 
than the Vink et al. predictions by 0.9 dex for HD 96715 
and 0.5 dex for HD 190429A. On this basis, Bouret et al. 
conclude that our understanding of O-star winds requires 



fundamental revision. The success here in reproducing their 
estimates shows that the required revision is simply the 
inclusion of a dynamically consistent transition to the ob- 
served photosphere. In addition, this theoretical confirma- 
tion of their low rates reinforces their conclusion that the 
evolutionary tracks of massive stars need to be recomputed. 

3.5. Mechanism 

The internal details of the RL models can be investigated 
in order to understand the sharp rise in J as one moves 
from late- to early- type O dwarfs. Most illuminating are 
the fractional contributions of various ions to gi in the layer 
containing the sonic point. Such calculations have been car- 
ried out for the models with logg — 4. These reveal that the 
high J's for T eff Z 40, 000K are due to the ability of high 
ions to absorb momentum from the emergent flux below the 
Lyman limit, as has long been understood (Castor, Abbott, 
& Klein 1976; Lamers & Morton 1976). Of dominant impor- 
tance is the changing ionization balance of iron, resulting 
in the effective driving ion Fe V replacing the ineffective ion 
Fe IV (cf. Vink et al. 2000). Thus, in the weak-wind domain 
at T eff = 32, 500K , Fe V lines contribute only 0.7% of g e 
at v = a, but this rises to 44% at 40,000K and reaches a 
maximum of 78% at 45,000K. 

These results show that the termination of strong winds 
in (T e ff, g)-space is determined by ionization balance in 
the reversing layer. In these models (L07a, Sect. 2. 5), this 
balance is matched at a representative point to that of the 
corresponding static TLUSTY atmosphere. 



4. Conclusion 

The aim of this paper has been to respond to an appar- 
ent problem for the theory of radiatively-driven winds aris- 
ing from the extraordinarily low $'s found for late- type 
Galactic O stars, as exemplified recently by M09. The prob- 
lem is lessened, though perhaps not entirely resolved, when 
the observed rates are compared to predictions made from 
first principles using the theory of dynamic RLs (LS70; 
L07a) rather than to the semi-empirical estimates of Vink 
et al. (2000). When this theory is used to define a weak- 
wind domain in (T e ff, g)-space, the M09 stars are found to 
be within or at the edge of this domain. It remains to be 
seen whether future developments in our understanding of 
low density supersonic winds will close the remaining gap 
between predicted and observed rates. 

A still outstanding issue is the huge disparity between 
the $y's and the observed estimates. It would be informa- 
tive if the Vink et al. code were used to fit the M09 stars 
individually and the emergent MC spectra compared with 
the observed spectra, as Abbott & Lucy (1985) did for £ 
Puppis. From experiments reported by M09, this consis- 
tency check is expected to fail, revealing that the $y's give 
rise to unacceptably strong wind lines. Such discrepancies 
would then provide the basis for revising the Vink et al as- 
sumptions about physical conditions at velocities ~ 0.5Uoo 
in these stars' winds. 

This substantial progress towards resolving the weak- 
wind problem and the success in reproducing the clumped- 
wind $'s of Bouret et al. (2005) strongly supports the sim- 
ple picture for the onset of winds presented in LS70: se- 
lective radiation pressure expels high atmospheric layers 



6 



Lucy: Mass fluxes 



causing a pressure inbalance in the photosphere resulting 
in an upwelling of matter from deeper layers. The natu- 
ral steady-state configuration is then the laminar transonic 
flow described by the theory of moving RL's (LS07; L07a), 
which merits still further development. 

Acknowledgements. I am grateful to J.S.Vink, A. dc Koter and 
T.L.Hoffmann for prompt answers to queries relevant to the refer- 
eeing of this paper. 
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